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ABSTRACT 


This  thesis  deals  with  a  simulation  of  a  missile  versus  target  engagement 
scenario.  After  deriving  simplified  transfer  functions  for  the  missile  seeker 
head,  missile  autopilot,  missile  dynamics,  and  tiirget  dynamics,  a  three 
dimensio’iial  simulation  is  developed  using  classical  proportional  navigation.  The 
scenario  is  simulated  using  state  variable  design.  A  forward  time  solution  of  the 
two  dimensional  problem  is  developed  which  is  converted  to  an  adjoint  model. 
The  adjoint  model  is  used  to  determine  the  optinr.al  time  to  initiate  simplified  and 
tactical  evasion  maneuvers  in  order  to  maximize  the  final  miss  distance. 


THESIS  DISCLAIMER 


The  reader  is  cautioned  that  computer  programs  developed  in  this 
research  may  not  have  been  exercised  for  all  cases  of  interest.  While  every 
effort  has  been  made,  within  the  time  available,  to  ensure  that  the  programs  are 
free  of  computational  and  logical  errors,  they  cannot  be  considered  validated. 
Any  application  of  these  programs  without  additional  verification  is  at  the  risk  of 
the  user. 
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1.  INTRODUCTION 


Since  the  early  years  of  the  Vietnam  War  the  surface-to-air  missile  has 
become  the  largest  threat  to  combat  aviators.  During  the  emergent  age  of  low 
intensity  conflicts,  a  goal  of  zero  percent  aircraft  losses  must  remain  a  high 
priority  to  strike  planners  and  decision  makers.  In  order  to  minimize  losses  to 
enemy  antiair  weapon  systems,  in  particular  surface-to-air  missiles,  these 
systems  must  be  generally  understood  and  tactics  necessary  to  defeat  these 
missiles  must  be  investigated. 

Aircraft  combat  survivability  is  deHned  as  "the  capability  of  an  aircraft  to 
avoid  and/or  withstand  a  man-made  hostile  environment."  [Ref.  1]  Hgure  1 
displays  the  basic  aircraft  survivability  relationships.  The  probability  of  an 
aircraft  being  hit  by  antiair  weapons  is  termed  Ph  which  is  referred  to  as  the 
susceptibility  of  the  aircraft  Susceptibility  is  calculated  by  the  product 

=  Pa  Pdit'Plgd 

where 

Ta  -  probability  the  threat  is  active  and  ready  to  engage  the 

aircraft 

Pnr  ~  l^'obability  the  aircraft  is  detected,  identified,  and  track^. 

Plgd  -  probability  the  threat  is  launched,  guided  and  either  hits  the 

aircraft  or  detonates  close  enough  to  cause  a  hit. 

The  vulnerability  of  an  aircraft  is  deHne^i  as  the  inability  to  withstand  the 

damage  caused  by  a  hostile  environment  It  can  be  measured  by  P]c/h  which  is 

the  probability  of  an  aircraft  kill  given  it  is  hit  by  hostile  fire.  This  leads  to  the 

relationship 

probability  of  loll suscepdbility'vulnerability,  (1.2) 
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Figure  1.  Aircraft  Survivability  Diagram 


or  in  ether  words, 

I^  =  I^Pk/h.  (1.3) 

It  is  obvious  then  that  the  probability  of  survival  is 

Ps»1-Pk.  [Ref.  1]  (1.4) 

Yherefoirt.  in  order  to  increase  the  survivability  of  an  aircraft,  the  probability 
of  kill  must  be  reduced.  One  way  that  this  can  be  achieved  is  by  decreasing 
Plcd.  The  probability  that  a  missile  is  successfully  guided  and  Jther  hits  an 
aircraft  or  detonates  close  enough  to  hit  it,  can  be  reduced  by  degrading  the 


guidance  and  control  system  of  the  missile.  The  employment  of  high 
acceleration  turns  and  random  maneuvers  are  examples  of  methods  that  can  be 
used 

This  research  develops  a  missile/target  simulation  progrnn  using 
classical  proportional  navigation.  A  tiiree  dimensional  model  is  produced 
assuming  a  dual  gimbal  axis  seeker  head.  Chapter  11  introduces  the  idea  of 
proportiobal  navigation  and  the  actual  guidance  law  is  developed.  In  Chapter 
m  the  transfer  functions  for  the  individual  missile  subsystems  are  determined. 
The  specifics  of  the  computer  simulation  geometry  are  discussed  leading  to  the 
actual  relationships  used  in  the  program.  Additionally,  a  two  dimensional 
forward  time  model  is  developdl  which  is  converted  to  an  adjoint  model.  The 
adjoint  model  is  a  useful  tool  for  analyzing  time  varying  systems.  It  is  used  in 
this  thesis  to  determine  optimal  miss  distance  parameters.  Chapter  W  consists 
of  actual  simulation  results.  Two  typical  scenarios  are  conducted  to  show  the 
capabilities  of  the  program  against  various  target  maneuvers.  The  adjoint 
model  is  then  used  to  determine  an  optimal  time  to  initiate  various  evasion 
maneuvers.  Conclusions  and  recommendations  follow  in  Chapter  V. 

All  computer  sirrulations  are  developed  and  conducted  using  the  Matrix 
Laboratory  (M/.TLAB)  and  the  three  dimensional  plots  are  generated  using 
die  Display  Integrated  Software  System  and  Plotting  Language  (DISSPLA). 


II.  BASIC  PROPORTIONAL  NAVIGATION 


A.  GENERAL 

Proportional  navigation  is  the  basic  guidance  law  used  in  the  majority  of  the 
threat  missiles  in  operational  use  today.  This  method  of  guidance  generates 
missile  acceleration  commands  proportional  to  the  line  of  sight  (LOS)  rate. 
Figure  2  shows  the  basic  parameters  and  geometry  associated  with  the  two 
dimensional  missile/target  engagement  problem  and  the  parameters  are  deflned 
below: 

Vt  =  target  velocity. 

Vm  =  missile  velocity. 

71  =  target  flight  path  angle. 

Ym  s  missile  flight  path  angle. 

R  s  missile  to  target  range. 

Rt  =  target  range. 

Rm  =  missile  range, 

a  =  line  of  sight  angle. 

0(  =  target  line  of  sight  angle. 

Cm  =  missile  line  of  sight  angle. 

B.  CONSTANT  BEARING  COURSE 

A  constant  bearing  course  is  one  where  the  line  of  sight  between  the 
target  and  the  missile  maintains  a  constant  orientation  in  space.  As  a  result, 
progressive  lines  of  sight  remain  parallel  to  each  other  as  the  engagement 
precedes.  Similarly,  the  line  of  sight  angle,  c,  remains  constant.  Hgure  3 


Figure  2.  Missile/Target  Geometry 


depicts  the  constant  bearing  course  idea.  As  long  as  there  is  a  positive  closing 
velocity  between  the  missile  and  the  target,  the  constant  bearing  course 
concept  will  ensure  an  intercept.  Proportional  navigation  uses  the  constant 
bearing  course  idea  by  driving  the  line  of  sight  rate,  d,  to  zero.  [Ref.  2] 


5 


C.  PROPORTIONAL  NAVIGATION  SCHEME 

Figure  4  depicts  the  basic  proportional  navigation  schme.  Assuming  that 
the  seeker  head  of  the  missile  follows  the  target,  the  transverse  acceleration 
perpendicular  to  the  lin'.  of  sight  will  equal  the  acceleration  of  die  R  vector  in 
that  direction.  Mathematically,  the  accele.'adon  of  R  is 

Ar  =!(R  +  <oxo)xR)iK+(2toxRxR+6xR)4  (2.1) 

where 

R  s  missileAarget  line  of  sight  vector. 


R  =  closing  rate  along  R. 

R  =  acceloation  along  R. 

(0  =  angular  rate  of  change  of  R  in  inertial 

space. 

Am  missile  acceleration  perpendicular  to  i.. 

At  =  target  acceleration  perpendicular  to  R. 
Ar  =  overdl  acceleration  of  R. 


At  this  point,  a  missile  acceleration.  Am*  equal  to  the  target  acceleration.  At,  will 
make  the  line  of  sight  parallel  to  its  original  direction.  As  long  as  R  remains 


Figure  4,  Vectorial  Proportional  Navigation  Scheme 
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along  R  ((0=0)  a  missile/target  impact  is  assured.  So,  the  transverse 
acceleration  conunand  is 

At-Am  =<bxR+2((0xR).  (2.2) 

Assuming  the  line  of  sight  rate  is  equal  to  the  angular  rate  of  change  of  R  in 
inertial  space,  equation  (2.2)  now  becomes 

At-A„  =  Rd+2Rd.  (2.3) 

D.  PROPORTIONAL  NAVIGATION  GUIDANCE  LAW 

In  the  classical  proportional  navigation  scheme,  the  missile  course  is  one  in 
which  the  rate  of  change  of  the  missile  heading  is  directly  proportional  to  the 
rate  of  rotation  of  the  line  of  sight  vector  from  the  missile  to  tlte  target.  As  a 
result,  this  course  change  is  intended  to  counteract  the  rot^on  of  the  line  of 
sight,  thus  returning  to  a  constant  bearing  course.  The  movement  of  the  missile 
and  target  cause  the  line  of  sight  to  rotate  resulting  in  a  differential  displacement 
between  the  missile  and  the  targe,  perpendicular  to  the  range  line.  Hgure  5 
depicts  this  geometry.  [Ref.  3]  The  proportional  navigation  guidance  law 
attempts  to  generate  an  acceleration  command,  Ac,  perpendicular  to  the  line  of 
sight 

Assume  a  gyro  stabilized  seeker  head.  If  there  is  no  torque  applied  to  the 
gyro,  the  seeker  will  not  rotate.  Assuming  the  seeker  tracks  the  target  the 
gyro  angle  will  follow  the  line  of  sight  Applying  the  equation  of  motion  for  a 
gyro  stabilized  seeker 

L  =  I(i)Q  a4) 

where 

L  =  applied  torque. 

DO  =  Spin  angular  velocity. 
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I  =  moment  of  inertia  of  the  gyroscope. 

£2  =:  rate  of  precession  of  the  gyroscope. 

Applying  this  to  the  case  when  the  seeker  head  tracks  the  target,  Q  is  then 
replaced  by  the  rate  at  which  the  gyro  is  torqued  in  space.  This  is  simply  6 
which  is  the  line  of  sight  rate.  [Ref.  4]  Thus,  equation  (2.4)  becomes 

Lsicoa.  (2.5) 


This  torque  is  in  turn  applied  to  the  control  surface  of  the  missile  leading  to 
the  relationship 

A„*kL  =  kIa)a  (2.6) 

where  k  is  a  constant  of  proportionally.  Referring  to  Figure  6,  a  relationship  is 
determined  for  Am  in  terms  of  the  rate  of  change  of  the  missile  flight  path  angle, 
Y..  Given  the  missile  velocity  vector  at  some  point  in  time,  Vn(t),  and  suppose 
the  missile  undergoes  an  xceleratitm,  A^i,  during  an  interval  of  time,  dL  The 
velocity  vector  is  then  displaced  and  is  represented  by  the  vector  Voi(t4dt). 
The  angle  the  vector  is  aaversed  is  simply  dYin>the  differential  missile  flight 
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figure  6.  Missile  Acceleration  Relationship 

path  angle.  [Ref.  4]  For  small  angles  (which  are  guaranteed  by  making  dt 
small)  the  following  relationship  is  obvious 

A,ndt  =  VB,dyn,.  (2.7) 

Dividing  equation  (2.7)  by  the  time  interval,  dt,  the  missile  acceleration  is 
defined  as 

C2.8) 

ComtK:  ,  ^nations  (2.6)  and  (2.8) 

V„t„  =  klaw.  (19) 

Dividing  through  by  Vm,  the  proportional  navigation  guidance  law  becomes 


or 

Ym=Na.  (2.11) 

Equation  (2.11)  represents  the  classical  proportional  navigation  equation  where 

Ym  =  rate  of  change  of  the  missile  heading. 

®  =  rate  of  change  of  the  line  of  sight. 

N  =  proportional  navigation  ratio. 

The  navigation  ratio  determines  the  sensitivity  of  the  missile  system.  A  high 
navigation  ratio  will  lead  to  rather  high  gains  resulting  in  large  missile 
commands  for  small  changes  in  the  line  of  sight  rate.  On  the  other  hand,  small 
values  for  N  will  lead  to  small  missile  commands  for  a  given  d.  Larger 
navigation  ratios  are  preferred  for  head  on  engagements  and  smaller  ones  are 
preferred  for  tail  chase  cases.  For  this  research  the  navigation  ratio  between 
three  and  five  was  chosen  [Ref.  4]. 
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in.  SIMULATION  DEVELOPMENT 


A.  INTRODUCTION 

The  development  of  the  actual  seeker  head  and  autopilot  transfer  functions 
is  discussed.  Ideal  missile  and  target  trajectory  equations  are  used  and  the 
state  equations  for  these  systems  are  generated.  Basic  geometric  relationships 
are  used  to  develop  the  missile  and  target  proportional  navigadcn  equations. 
Once  these  equations  are  determined  for  the  two  dimensional  problem,  the  state 
equations  are  then  augmented  to  a  three  dimensional  problem.  The  continuous 
state  equations  for  the  seeker  head,  autopilot,  missile  kinematics,  and  target 
kinematics  are  cpnverted  to  equations  appropriate  for  digital  simulation. 
Hnally,  the  two  dimensional  forward  time  model  is  developed  and  converted  to 
an  adjoint  model. 

B.  MISSILE  SUBSYSTEM  DEVELOPMENT 

1.  Subsystem  Layout 

A  basic  functional  block  diagram  of  a  generic  tactical  missile  system  is 
shown  in  Figure  7.  Although  the  exact  configuration  and  description  of  each 
component  depends  on  many  different  factors,  an  attempt  will  be  made  to 
develop  a  basic  system  sufficient  for  simplified  simulation.  Figure  8  shows  the 
typical  physical  location  of  the  various  missile  subsystems. 

2.  Seeker  Head  Development 

The  seeker  head  can  simply  be  thought  of  as  the  eye  of  the  missile.  It 
is  able  to  detect,  acquire,  and  track  a  target  by  sensing  some  unique 
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Figure  7.  Missile  Subsystem  Block  Diagram 
characteristic  of  the  target  itself.  This  usually  consists  of  the  radiation  or 
reflection  uf  energy  by  the  taiget 

A  seeker  with  a  namw  field  of  view  will  be  used.  Hgure  9  shows  a 
basic  gimballed  seeker  head  configuration.  Here,  the  actual  seeker  is  mounted 
on  a  gimballed  platform  and  it  maintains  the  target  within  the  field  of  view  by 
rotating  the  platfontL  The  inertial  rotation  rate  of  the  line  of  sight  provides  the 
missile  with  the  required  tracking  information.  [Ref.  S] 

Hgure  10  displays  a  diagram  for  an  actual  seeker  head  where 
P  3  seeker  head  gimbal  angle. 

The  control  torque  to  the  seeker  results  in  the  following  equation  of  motion 


T»Ip 

where 

T 

s  applied  torque. 

I 

s  moment  of  inertia. 

P 

-  angular  acceleration. 
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Figure  8.  Basic  Missile  Layout 

Solving  for  p 

P-I.-k,(B-<,)-k2p  (3.2) 

« -k2p-kiP+k|*7.  (3.3) 

The  constants  ki  and  k2  are  determined  by  the  time  constant  of  the  seeker  head 
system. 


Now  taking  the  Laplace  Transform  of  equation  (3.3)  and  assuming 


£{p}.£{-kjP-k,p+k,o}  (3.4) 
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zero  initial  conditions 


s^(s)s-k2Sp(s)-kiP(s)+kia(s)  (3.5) 

8^(s)+k2Sp(s)+kip(s)  =  kiO(s).  (3.6) 

Finally,  the  tran.sfer  functiOT*  beco<nes 


M. 


EL 


EL 


(3.7) 


o(i)  s^  +  k2i+k|  (s+l/Tjhf  * 

For  this  simulation  the  time  constant  of  the  seeker  head,  t,|,,  was 
chosen  to  be  0.1  second  which  is  a  good  approximation  of  a  real  world  system. 
From  equation  (3.7)  the  constants  ki  and  kz  are  defined  in  terms  of  this  time 
constant  where 


-100,  (3.8) 

k2*2(l/T,„)-20.  (3.9) 
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Figure  1 1  depicts  the  signal  flow  graph  of  the  seeker  head  system. 
From  this  diagram  the  continuous  time  state  equation  of  the  form 

Xsh  =  Axsh+Bush  (3-10) 

is  easily  determined.  Selecting  the  state  vector  to  be 


and  the  input  =  a,  equation  (3.10)  becomes 


(3.11) 


'su 


=  [-i°00  -2oh*\°i. 


“sh- 


(3.12) 


Figure  11.  Seeker  Head  Signal  Flow  Graph 
3.  Guidance  and  Autopilot  Development 

The  guidance  law  within  the  ntissile  system  determines  the  best 
trajectory  for  the  missile  based  on  the  missile  position,  target  position,  missile 
capability  and  the  desired  objectives.  A  command  is  sent  to  the  autopilot  which 
determines  the  control  (i.e.,  actuator  position  and  thrust)  necessary  to  perform 
such  a  command.  In  the  proportional  navigation  guidance  law,  the  missile 
commando  are  generated  in  order  to  change  the  missile  flight  path  rate 
proportionally  to  the  missile  to  target  line  of  sight  rate. 
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In  this  simulation  a  very  simplified  autopilot  will  be  used.  Referring  to 
Figure  12,  the  transfer  function  is  easily  developed.  The  angle  of  attack  of  the 
missile  is  assumed  to  be  zero,  so  the  velocity  zf  the  missile  will  be  aligned  with 
the  missile  center  line  at  an  angle,  Ym*  Applying  a  torque  to  the  missile  about 
the  center  of  gravity  results  in  tl  following  equation  of  motion 

TMm  ~  ^cgY  (3.13) 

where 

=  commanded  torque. 

leg  =  moment  of  inertia  about  the 

center  of  gravity. 

Y  =  angular  acceleration  of  the 

missile  flight  path. 


Figure  12.  Basic  Missile  Layout 
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Solving  for  Y, 


y  =  ^  -  +  kNP  (3.14) 

leg 

where  k  is  determined  by  the  slowest  time  constant  of  the  missile/autopUot. 

Now  taking  the  Laplace  Transform  of  equation  (3.14)  and  assuming 
zero  initial  conditions 

•£•{?»,}  =  i-{fcrm+kNp}  (3.15) 

s*Tfm(s)  =  -ltSTm(s)+l£Nsr„(s)  (3.16) 

(s^  +  ks)Y„(s)  =  lcNsp(s).  (3.17) 


Hnally,  the  transfer  function  becomes 


(3.18) 


Tp(0-  fcN 
P(s)  s*k 

where  k  is  defined  above.  The  autopilot  time  constant  ,t^,was  chosen  to  be  1.0 
second,  so 


k  =  l/T,p*1.0. 


(3.19) 


Figure  13  shows  the  signal  flow  graph  for  the  missile  autopilot  system. 
From  this  diagram  the  continuous  state  equation  is  easily  determined.  Selecting 
the  state  vector  to  be 

(3^) 

and  the  input  UqpN^,  the  state  equation  becomes 

*«p*[-*l]*ap+[l]«ap-  (3^1) 
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4.  Missile  and  Target  Dynamics 


The  signal  flov^  graphs  for  the  target  and  missile  dynamics  are  shown 
in  Rgure  14  where  A„=Vn,Y,n.  At  is  selected  independently.  The  state 
vectors  are 


and 


*tn  “ 


ym 

ym 


(3.22) 


Xt 

X»  =. 

yt 
Lyt 


(323) 


With  and  ut=[““ 

L“inyJ  [“ty. 


,  the  state  equations  become 
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and 
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*t  = 


o 

o 

o 

o 

o 

0  0  0  0 

Xt  + 

1  0 

0  0  0  1 

1 

0  0 

0  0  0  0. 

0  1 

u. 


(325) 


5.  Overall  System  Layout 

Hgure  15  shows  the  total  system  signal  flow  graph.  The  seeker  head 
angle  rate,  is  the  estimate  of  the  line  of  sight  rate,  dr. 

6.  Three  Dimensional  State  Definitions 

Assuming  that  the  missile  is  roll  stabilized  and  controllable  in  both  pitch 
and  yaw,  the  state  equations  for  the  missile  subsystems  are  extended. 
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The  seeker  head  state  matrix  now  becomes 


Ppiich 


Xsh  = 


026) 


Kyaw 


and  the  continuous  state  equation  is  obtained  by  simply  augmend'ij  the  two 
dimensional  problem. 

The  resulting  equation  becomes 


Xjii  ss 


■  0  1 

-100  -20 
0  0 


0  0  -100  -20 J 


100  0 
0  0 
0  100 


where  u^i, 


The  autopilot  state  matrix  now  becomes 


[Yni.pitch1 

tra_y«w  J 


(327) 


(328) 


and  the  continuous  state  equation  is  again  simply  obtained  from  an 
augmentation  of  the  two  dimensional  equation.  The  resulting  equ^on  for  the 
autofrilot  is 


(329) 
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where 


LHyawJ 

For  the  missile  and  target  dynamics  equations  the  state  equations  are 
simply  augmented  with  vertical  position  and  vertical  velocity  states. 

These  equtttions  are 


X  = 


and  the  state  variable  dynamics  become 


(330) 
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(331) 


where  u 


L“xj 


,  the  force  per  unit  mass  acting  on  the  missile  or  target. 
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C.  MISSILE  AND  TARGET  GEOMETRY 
1.  Velocity  Relationships 

Hgure  16  shows  the  vectorial  relationship  of  the  missile  and  target 
velocity.  From  this  figure 

and 

V,  =  t/V+V^^-  (3-33) 


2.  Line  of  Sight  Angled 

In  order  to  simulate  the  scenario  in  three  dimensional  space  a  line  of 
sight  angle  in  pitch  and  yaw  must  be  defined.  Figure  17  depicts  these  angular 
relationships.  Given  the  target  and  missile  position  in  Cartesian  coordinates  the 
line  of  sight  angles  are  defined  as 

®piich  *  j  (3*34) 
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Figure  17.  Pitch  and  Yaw  Line  of  Sight  Angle  Definitions 
and 

Oy^»tan“‘[(y,-y„)/(xt-x„)l.  (335) 

Figure  18  shows  the  target  and  missile  line  of  sight  ang'e  definitions  and  these 
angles  are  expressed  mathematically  as 


®in_pildi  “  ^  [*in  !  ]♦ 

(336) 

®ni_yiw  ”  ^ 

(337) 

and 

'  ' 

(338) 

Oty^-tan-^lyt/Xt]. 

(339) 
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Figure  18.  Target  and  Missile  Line  of  Sight  Angle  Definitions 
3.  Flight  Path  Angles 


The  missile  and  target  flight  path  angles  are  also  defined  in  three 
dimensions.  Figure  19  depicts  the  geometric  definitions.  Mathematically, 


Figure  19.  Target  and  Missile  Flight  Path  Angle  Definitions 
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(3.40) 


Ym.pilch  ~  [ymz  !  Yny^  j» 

Yin_y«w~*^  [Yny/^inx]»  (3.41) 

and 

Yt_pteh  =  tan-'  [  ]•  (3-42) 

Yi.y«,  =  tan-'[V^/V^].  (3.43) 

4.  Velocity  and  Acceleration  in  the  Pitch  Plane 

The  missile  is  controlled  in  three  dimensional  space  by  generating 
acceleration  commands  in  two  orthogonal  planes.  These  planes  are  defined  as 
the  pitch  plane  and  the  yaw  plane.  The  yaw  plane  is  taken  as  simply  the 
horizontal  XY  plane  and  the  pitch  plane  is  the  orthogonal  plane  rotated  by  the 
angle  Oyaw.  Figure  20  depicts  the  geometric  definition  of  the  pitch  plane.  Given 
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the  missile  velocity,  and  the  angles  Cyaw  and  Ym_yaw«  the  velocity  in  the  pitch 
plane  is 

^in_iriich  “  — <Tyj^  (3.44) 

From  the  basic  relationship  for  the  missile  acceleration  the  missile 
pitch  acceleration  is  defined  as 

An,_pitdi  ~  Yn.piich  "‘tm.pildi*  (3.45) 

This  acceleration  vector  is  then  broken  down  into  Cartesian  coordinate  system 
components.  Assuming  that  the  pitch  acceleration  vector  is  perpendicular  to 
the  vertical  line  of  sight  vector  between  the  missile  and  target,  Hgure  21  shows 
the  orientation  of  the  acceleration  components.  From  this  Hgure  the  following 
relationships  are  obtained 

^pilch  ) ^ytw*  (3.46) 

yin_piicfa  *  ”{ Ain_pildi ' ^pitch  ) ^yiw»  (3.47) 


Figure  21.  Pitch  Plsne  Acceleration  Components 


aad 

^m_pUch  ~  ■^in_pitch  (3.48) 

5.  Velocity  and  ^acceleration  in  the  Yaw  Plane 

The  missile  yaw  plane  is  depicted  in  Figure  22.  Given  the  missile 
velocity  and  the  angle  Tm_puch>  the  velocity  in  the  yaw  plane  is  singly 

^n_yaw  “  V|n ‘C®*Yin_piich*  (3.49) 


Again,  from  the  basic  missile  acceleradoc  relationship,  the  missile  yaw 
acceleration  is  defined  as 

*  YB_3r«w’Yai^3ftw  (3.50) 

The  missile  yaw  acceleration  components  are  depicted  in  Figure  23  and  are 
given  mathematically  as 

*in_3r«w  “ (3.51) 
and 
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(3^2) 


yin_y«w  =  cbsCTya^. 


Figure  23.  Yaw  Plane  Acceleration  Components 
6.  Total  Missile  Acceleration  Components 

The  input  to  the  missile  state  equation  is  the  vector  consisting  of  the 
three  Cartesian  components  of  the  overall  missile  acceleration  vector.  Given 
the  missile  acceleration  conqwnents  in  both  the  pitch  and  yaw  planes,  the  total 
missile  acceleration  components  are 


(333) 

y«  “  yin_pildi  ■*"ym_ytw» 

(334) 

and 

(335) 

Hnally,  the  total  missile  acceleration  is 

(336) 

7.  Target  Acceleration  Components 

In  this  simuiatibn  the  target  acceleration  conr.ponents  consist  of  the 

three  Cartesian  components  of  the  overall  target  acceleration  vector.  These 
components  are  used  as  the  input  to  the  target  idhematic  equations.  The  ininit 


vector  is  determined  by  the  type  of  target  maneuver  desired.  The  overall  target 
acceleration  is  simply 

At  *  V^t^  +  yt^+zt^.  (3.57) 

8.  Closing  Velocity  and  Time  To  Go 

From  equation  (3.48)  the  velocity  of  the  target  in  the  pitch  plane  is 

sinqjly 

yt_pitch  “  Yt  “^ysw)*  (338) 

Now,  the  range  rate,  R,  is  found  by  simply  projecting  the  missile  and  target 
pitch  plane  velocities  along  the  vertical  line  of  sight  from  the  missile  to  the 
target  Mathematically 

R  =  •COs(Yl_pilgli  "'<ypitch)*“  •COs(Yni_pilch  “®piJdi)*  (339) 

Knowing  that  the  range  rate  is  the  negative  of  the  closing  velocity,  the  time  to 
go  is  simply 

Tmie-to-go  =  -— ^ —  =  -^,  (3.60) 

Yclosing 

where  R  is  the  range  from  the  missile  to  the  target. 

D.  DIGITAL  SIMULATION  USING  STATE  SPACE  METHOD 
1.  Discrete  State  Equation  Definition 
Given  the  continuous  state  equations 

x(t)  =  Ax(t)+Bu(t),  (3.61) 

y(t)-Cx(t)+Du(t),  (3.62) 

a  more  convenient  form  for  digital  simulation  is 

x(k  +  l)*Ox(k)+ru(k),  (3.63) 

y(k+l)  =  Cx(k)+Du(k).  (3.64) 


The  discrete  state  matrices  are  defined  as 

O  =  =  I  +  A  •  DT+ (1  /  2!)A^  •  DT^  +  (1  /3 !)  A^  •  DT^ +.. .,  (3.65) 

r  =  [lDT+(l/2!)ADT^+(l/3!)A2DT^+...]B,  (3.66) 

and  for  u(t)  constant  over  [(k+l)-k]  DT. 

2.  Missile  Subsystem  Discrete  State  Equations 

The  missile/target  scenario  is  simulated  digitally  using  the  MATLAB 
software  package.  A  built  in  system  function  is  used  to  convert  the  continuous 
state  equations  to  discrete  time  state  equations.  The  discrete  state  equations 
used  in  the  simulation  are  defined  as  follows  for  the  seeker  head,  autopilot, 
missile  dynamics  and  target  dynamics  respectively; 

Xjh  (k  + 1)  -  <I>sh5tsh  (*c)+ rsh“sh  (k). 

and 

Xt(k  +  l)=:OiXt(k)+rtUt(k). 

E.  FORWARD  TIME  AND  ADJOINT  MODEL 

To  accurately  and  simply  assess  miss  distances  for  various  target 
maneuvers,  a  two  dimensional  adjoint  model  is  used.  Brst,  a  forward  time 
model  is  developed  for  the  two  dimensional  missile/target  engagement  Figure 
24  shows  the  block  diagram  of  the  forward  time  model.  All  system  inputs  are 
converted  to  impulses  for  subsequent  conversion  to  the  adjoint  model. 

The  forward  time  mode  is  used  to  show  the  miss  distance  due  U)  heading 
error  (where  heading  error  is  defined  as  the  error  in  missile  heading  from  a 


(3.67) 

(3.68) 

(3.69) 

(3.70) 
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collision  course)  and  the  miss  distance  due  to  target  maneuver.  The  missile  to 
target  range  is  defined  as 

R  =  Vctime-to-go  (3.71) 


where 


time“to-go  =  Tfinai-t.  ,  (3.72) 

Tlds  model  shows  the  miss  distance  throughout  the  course  of  the  engagement . 
The  last  value  of  miss  distance  is  the  miss  distance  at  Tfin«i. 

An  adjoint  model  is  used  to  analyze  linear  time  varying  systems  [Ref.  6]. 
The  major  advantage  of  the  adjoint  technique  is  that  the  miss  distance  for  all 
final  times  is  generated  in  one  run  rather  than  the  multiple  runs  required  when 
using  a  Monte  Carlo  simulation.  Following  the  rules  for  adjoint  construction 
[Ref.  6],  Figure  24  is  converted  to  the  adjoint  model  depicted  in  Hgure  25. 


Figure  24.  Forward  Time  Model 
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Figure  25.  Adjoint  Model 


IV.  SIMULATION  RESULTS 

A.  OVERVIEW 

This  section  presents  the  results  of  the  computer  simulations.  The  scenario 
is  conducted  for  several  different  initial  conditions  and  with  several  different 
target  maneuvers.  The  following  assumptions  are  made  throughout: 

1.  The  missile  is  limited  to  20  g's  in  either  plane. 

2.  The  acceleration  due  to  gravity  is  ignored. 

3.  The  target  is  Capable  of  instantaneous  acceleration. 

4.  The  target  is  limited  to  a  maximum  of  8.0  g's. 

5.  The  maximum  missile  speed  is  3000  feet  per  second. 

6.  The  maximum  target  speed  is  ISOO  feet  per  second. 

7.  The  proportional  navigation  constant  is  4. 

8.  The  missile  is  pointed  toward  the  target  at  launch. 

B.  THREE  DIMENSIONAL  ENGAGEMENTS  AGAINST 
VARIOUS  TARGET  MANEUVERS 

1.  Scenario  1:  A  Constant  Velocity  Target 

The  simulation  is  initially  conducted  with  the  missile  engaging  a 

constant  velocity  target  Figure  26  depicts  the  scenario  geometry.  The  initial 


conditions  are 

xm(0) 

s?  Ofeet 

ym(0) 

~  0  feet 

2m(0) 

®  Ofeet 

vmx(0) 

ss  3000  fps 
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Figure  26.  Constant  Velocity  Target  Scenario  Geometry 

vmy(O)  =  Ofps 
vmz(0)  =  0  fps 
xt(0)  =  30000  feet 

yt(0)  =0fect 
zt(0)  =  500  feet 

vtx(0)  =  0  ^ 
vty(0)  =  1000  Ips 
vtz(0)  =  Ofps. 

Hgure  27  depicts  the  scenario  that  reflects  these  initial  conditions. 
In  this  case,  the  aircraft  is  attempting  to  stay  on  the  outer  edge  of  the  missile 
envelope. 
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Figure  27.  Tactical  Scenario  for  Constant  Velocity  Target 
2.  Scenario  2:  A  Constant  Acceleration  Target 

la  this  scenario,  the  target  is  accelerates  at  a  constant  rate  of  15 


feet/sec/sec  in  the  y  direction.  Figure  28  shows  the  scenario  geometry.  The 


initial  conditions  arc 

xm(b) 

-  Ofeet 

ym(0) 

-  21000  feet 

zni(C) 

s  Ofeet 

vmx(0) 

»  2100fps 

vmy(0) 

=  -2100  fps 

vmz(0) 

*  O  fps 

xt(0) 

-  21000  feet 

>t(0) 

=  Ofeet 

«(0) 

500  feet 
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Figure  28.  Constant  Acceleration  Target  Scenario  Geometry 


vtx(O)  =  Ofps 
vty(0)  =  1000  ^ 
vtz(0)  =  Ofps. 
xddt(0)  =  0fps2 
yddKO)  =  15fps2 
zddt(0)-0f^2 

Pd(0)  *  -4r. 

Figure  29  depicts  the  scenario  that  reflects  these  icitial  conditions.  In 
this  case  there  is  overlapping  missile  envelope  coverage  and  the  aircraft  is 
attempting  to  strike  a  target  located  within  this  region.  The  missile  is  launched 
when  the  target  enters  the  SNM  radius. 
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Figure  29.  Tactical  Scenario  for  Constant  Acceleration  Target 
3.  Scenario  3:  A  Two  Dimensional  Target  Acceleration 


Hgures  26  and  27  depict  the  initial  conditions  for  this  scenario.  The 


target  performs  a  6.5  g  level  turn  toward  the  missile  at  a  range  of  12^000  feet 

The  target  acceleration  is 

Xj  «* -6. 5  •  32. 2  •  sin^Yi  y«w ) 

(4.1) 

y.  =6.5*j2.2*cos^Yt^y*w) 

(4.2) 

Z^ssO.O. 

(4.3) 
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4.  Scenario  4:  A  Three  Dimensional  Target  Acceleration 

Figures  28  ard  29  depict  the  initial  conditions  for  this  scenario.  The 
target  performs  a  three  dimensional  maneuver  toward  the  missile  and  into  the 
vertical  at  a  range  of  12,000  feet  The  target  acceleration  is 


“  “6.5 *32.2 ‘Sin^Yi^y^j 

(4.4) 

yj  =  6.5'32.2*cos^Yi_yg^  j 

(4.5) 

~ 6.5‘32.2*cos^Yi  pitch)* 

(4.6) 

5.  Three  Dimensional  Simulation  Results 

Hgures  30-75  display  the  results  of  the  three  dimensional  simulations. 
Hgures  30-40  relate  to  the  constant  velocity  target  scenario.  Figures  41-52 
relate  to  the  constant  acceleration  scenario.  Hgures  53-64  relate  to  the  two 
dimensional  target  acceleration  scenario.  Hgures  65-75  relate  to  the  three 
dimensional  target  acceleration  scenario. 


FEET/SEC''2 


0  2  4  6  8  10  12 

TIME 


Figure  31.  Scenario  1:  Missile  Acceleration  jAnl 
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Figure  32.  Scenario  1:  Missile  Pitch  Flight  Path  Angle  Ta^pitch 


Figure  33^  Scenario  1:  Missile  Pitch  Flight  Path  Angle  Rate  7Bj)itch 
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r^^EKER  HEAD  PITCH  ANGLE  VS  TIME 


Figure  37.  Scenario  1:  Seeker  Head  Pitch  Angle  Rate  PpKch 
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Figure  38.  Scenario  1:  Seeker  Head  Yaw  Angle  P, 


Figure  39.  Scenario  It  Seeker  Head  Yaw  Angle  Rate  p. 
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Figure  44.  Scenario  2:  Missile  Pitch  Flight  Path  Angle  pjich 
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Figure  45.  Scenario  2:  Missile  Pitch  Flight  Path  Angle  Rate  Ym^itcb 
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Figure  48.  Scenario  2:  Seeker  Head  Pitch  Angle  Ppit^h 
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Figure  51.  Scenario  2:  Seeker  Head  Yaw  Angle  Rate 
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Figure  53.  Scenario  3:  Missile  to  Target  Range 


Figure  54.  Scenario  3:  Missile  Acceleration  \A^\ 
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Figure  55.  Scenario  3:  Target  Acceleration  |A|| 


Figure  56.  Scenario  3:  Missile  Pitch  Flight  Path  Angle  7^  pitch 


RADIANS  M  RADIANS/SEC 


RADIANS/SEC 


TIME 


Figure  60.  Scenario  3:  Seeker  Head  Pitch  Angie  ppn^ii 
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Figure  64.  Scenario  3:  Three  Dimensional  Plot 
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Figure  68.  Scenario  4:  Missile  Pitch  Flight  Path  Angle  Rate  tmjjitch 
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Figure  70.  Scenario  4:  Missile  Yaw  Flight  Path  Angle  Rate  Tm^w 
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Figure  71.  Scenario  4:  Seeker  Head  Pitch  Angle  Ppucii 


Figure  72.  Scenario  4:  Seeker  Head  Pitch  Angle  Rate  Ppit^h 


2'400.0 


C.  MISS  DISTANCE  ASSESSMENT 

A  two  dimensional  adjoint  model  is  used  to  investigate  the  missile  to  target 
miss  distance  due  to  a  variety  of  evasion  maneuvers.  The  optimal  time  to 
initiate  these  maneuvers  is  determined. 

Three  maneuvers  are  used: 

1.  Step  Target  Acceleration. 

2.  Barrel  Roll. 

3.  SpUt'S*. 


Additional  variables  are 

closing  velocity. 

= 

target  acceleration  magnitude. 

V„HE 

s 

missile  velocity  due  to  heading 

exTor. 

Wj 

* 

maneuver  frequency. 

L 

* 

half  period  of  Split  'S’. 

The  following  initial  conditions  are  constant  throughout  all  simulations: 

N 

s 

4X). 

Tf 

= 

6.0  seconds. 

Vc 

ss 

2500  feet  per  second. 

VnflE 

ss 

0.0  feet  per  second. 

3 

1.0, second. 

3 

0.1  second. 

1.  Step  Target  Acceleration 

The  step  acceleration  is  simulated  by  multiplying  a  step  input  signal  by 
the  magnitude  of  the  target  acceleration,  iif  Target  acceleration.'^  of  1.0,  6.0, 
and  8.0  g's  are  simulated. 
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2.  Barrel  Roll 

A  two  dimensional  Barrel  Roll  maneuver  can  be  represented  by  a 
sinusoid  of  frequency  cot.  A  shaping  Hlter  is  used  to  simulate  the  maneuver. 
[Ref.  7]  Rgure  76  depicts  a  typical  Barrel  Roll  maneuver  and  Figure  77  shows 
the  shaping  filter  equivalent  where  (0i  equals  1.0  radian/second.  The  filter  is 
converted  to  a  block  diagram  which  is  simulated  using  state  variable  design. 
The  block  diagram  is  shown  in  Hgure  78.  A  4.0  and  6.0  g  Barrel  Roll  maneuver 
are  simulated. 


B  ARREL  ROLL 

EASY  360*  ROLL 


Figure  76.  Typical  Barrel  Roll  Maneuver 


Uc 


1/(0 

t 

1  +  S  ^  /  CO  2 

t 

•• 

yt 


u.,  »  White  Noise  With  Power 
Spectral  Density  q.(0 
q^t)=1,,2/Tf  OstsTf 

"  0, Otherwise 


Figure  77.  Shaping  Filter  Equivalent  of  a  Barrel  Roll 
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1 
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S 

Figure  78.  Barrel  Roll  Block  Diagram 

3.  Split  ’S’ 

A  two  dimension;J  Split  maneuver  is  represented  by  a  periodic 
square  wave  with  period  2L.  A  shaping  filler  is  also  used  to  simulate  the 
maneuver.  [Ref.  7]  Hgure  79  depicts  a  typical  Split  'S'  maneuver  and  Figure 
80  shows  the  shaping  filter  equivalent  where  L  equals  1.0  second.  The  block 
diagram  of  this  filter  is  shown  in  Hgure  81.  A  4.0  and  6.0  g  Split  'S'  maneuver 
are  simulated. 
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4.  Forward  Time  and  Adjoint  Model  Simulation  Results 

Figures  82-95  show  the  results  of  the  forward  time  and  adjoint  model 
simulations.  Figures  82-87  show  the  results  of  the  step  acceleration  maneuver. 
Figures  88-91  show  the  results  of  the  Barrel  Roll  evasion  manetwer.  Figures 
92-95  show  the  results  of  the  Split  'S’  evasion  maneuver. 


Figure  79.  Typical  Split  'S'  Maneuver 


Us 


Uj  -  White  Noise  With  Power 
Spectral  Density  q.(t) 

0^‘STf 

•  0  .otherwise 


Figure  80.  Shaping  Filter  Equivalent  of  Split  ’S*  Maneuver 


Figure  81.  Split  'S*  Block  Diagram 


Tune 


Figure  85.  Adjoint,  6  g  Step  Acceleration 
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Figure  86.  Forward  Time,  8  g  Step  Acceleration 


Figure  87.  Adjoint,  8  g  Step  Acceleration 
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Time 


Figure  88.  Forward  Time,  4  g  Barrel  Roll 


Figure  89.  Adjoint,  4  g  Barrel  Roll 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 
A.  CONCLUSIONS 

Qeariy,  the  classical  proportional  navigation  model  of  the  missile  guidance 
and  control  system  is  sufficient  for  simple  trajectory  study.  The  use  of  the 
seeko*  head  angle  rate,  as  the  estimate  of  the  actual  line  of  sight  rate,  d, 
gives  accurate  enough  information  for  target  tracking  in  a  benign  environment. 
The  three  dimensional  model  of  the  missile  tracks  a  maneuvering  target  in  one, 
two,  or  three  dimensions  with  enough  accuracy  to  cause  an  aircraft  kill.  A 
navigation  ratio  of  4.0  generates  sufficient  missile  acceleration  commands  for 
the  sinq)lified  scenarios  studied. 

The  adjoint  model  proves  to  be  a  very  useful  indicator  of  miss  distance 
properties.  With  the  step  acceleration  evasion  maneuver,  the  optimal  time  to 
initiate  it  is  about  2.0  seconds,  or  5000  feet  (closing  velocity  is  2500  feet  per 
second),  prior  to  missile  impact  The  miss  distance  increases  as  the  target’s  g 
loading  increases.  The  maximum  miss  distance  of  140  feet  is  achieved  with  an 
8.0  g  maneuver. 

l^th  the  Barrel  Roll  maneuver,  the  optimal  time  to  initiate  it  is  about  3.5 
seconds,  or  8750  feet,  prior  to  missile  impact  The  miss  distance  again 
increases  with  increased  target  g  loading.  The  maximum  miss  distance  of  50(X} 
feet  is  achieved  with  a  6.0  g  maneuver. 

With  the  Split  'S'  maneuver,  the  optimal  time  to  initiate  it  is  about  1.8 
seconds,  or  4500  feet  prior  to  missile  impact  The  miss  distance  increases  with 
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increased  target  g  loading.  The  maximum  miss  distance  of  2600  feet  is 
achieved  with  a  6.0  g  maneuver. 

Geariy;  the  maximum  miss  distance  and  optimal  time  to  maneuver  depends 
on  the  evasion  maneuver  used.  In  general,  it  is  recommended  that  a  maximum 
g  evasion  maneuver  be  initiated  when  the  missile  is  between  1.0  and  1.5  miles 
prior  to  impact 

B.  RECOMMENDATIONS  FOR  FURTHER  STUDY 

The  simulations- developed  in  this  thesis  provide  a  very  good  generic  model 
which  can  be  easily  modified  for  more  complex  missile  systems.  The  programs 
are  easily  adjusted  for  systems  with  different  system  time  constants  and 
navigadon  ratios. 

The  problems  associated  with  measurement  noise  and  seeker  head 
measurement  errors  might  be  investigated  to  see  their  role  in  miss  distance 
evaluation. 

It  is  also  recommended  that  target  Electronic  Counter  Measures  (ECM)  be 
added  to  the  simulation  of  the  evasion  maneuvers  to  see  their  effects  on  miss 
distance. 


APPENDIX  A-THREE  DIMENSIONAL  PROGRAM 

dear 

clg 

% 

%  LT  Frands  C.  LukenbiU,  USN 

%  08  August  1990 

%  This  program  simulates  a  3  dimensional  target/missile 
%  engagement  scenario  using  classical  proportional  navigation. 

%  . 

%  DEHNE  STATES 

% 

%  Missile 
% 

%  ms=[xm  •  missile  x  coordinate 

%  xdm  -  missile  velodty-x  direction 

%  ym  •  missile  y  coordinate 

%  ydm  •  missile  velocity-y  direction 

%  zm  '  missile  z  coordinate 

%  zdm  missile  velocity-z  direction] 

% 

AM=[0  10  0  0  0 
000000 
000100 
000000 
oodool 
000000]; 

BM=[000 

100 

000 

010 
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000 
0  01]; 


Seeker  Head 

beta=[beta  pitch  -  seeker  head  pitch  angle 
betad  pitch  -  seeker  head  pitch  angle  rate 
beta  yaw  -  seeker  head  yaw  angle 
betad  yaw  •  seeker  head  yaw  angle  rate] 


AS=[  0  10  0 
-100  -20  0  0 
0  0  0  1 
0  0-100-20]; 

BS=[  0  0 
100  0 
0  0 
0  100]; 


% 

%  Aut 
% 

%  gar 
% 

* 

A?=l-1  0 
0-1]; 

BP=ri  0 
0  1]; 


Autopilot 

gammadm=[gammadni  pitch  -  pitch  body  angle  rate 
gammadm  yaw  -  yaw  body  angle  rate  ] 


Target 


t«s[xt  -  target  x  coordinate 


84 


%  xdt  -  target  velocity-x  direction 

%  yt  -  target  y  coordinate 

%  ydt  -  target  velocity-y  direction 

%  zt  >  target  z  coordinate 

%  zdt  -  target  velocity-z  direction] 

% 

AT=[0  10000 
000000 
000100 
000000 
000001 
000000]; 

BT-[000 
100 
000 
010 
000 
00  1]; 

% 

%  DISCRETE  REPRESENTATION 

% 

dt=.01; 

[phis.dels]-  c2d(AS3S,dt); 

[phinMlelinjs  c2d(AM3M,dt); 

[phia3ela]=  c2d(AP3P,dt); 

[phit,delt]=  c2d(AT3T,dO; 

tfinal-  15.0; 
kinax=tEna]/dt -i-l; 

% 

%  INITIALIZE  VARIABLES 

% 
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%  navigation  ratio 
% 

nr=54.0; 

NR=[nr  0 
0  nr]; 

% 

%  seeker  head  angles  and  angle  rates 

% 

beta^tchO  =0.0; 
bctad4>itch0=0.0; 
beta_^wO  =0.0; 
betad_yawO  =0.0; 

beta(:,l  )=P>cta_pitchO 
betad^tchO 
bela^awO 
betad^avO  I; 

% 

%  initial  missile  body  angle  rates 
% 

gsiininadmLjntchOsO.O; 
gaimmadm^awO  *0,0; 

g:  iminadm(:,l  )=[  gaminadmjntchO 

}  ganunadnLO'awO  ]; 


%  initial  missile  states 


xin0=0.0; 

x|ini0=:3000.0; 

ymO  =0.0;; 

ydmOO.0; 

zniD=O.0; 

zldinO=0.0; 
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ms(:,lM^unO 

xdmO 

ymO 

ydmO 

zmO 

zdmO]; 

% 

%  initial  target  states 
% 

xtO*  30000.0; 
xdi0=0.0; 
ytO  =0.0; 
ydtO=  1000.0; 
ztO  =500.0; 
zdt0=<).0; 

ts(:,l)=[xtO 

xdtO 

ytO 

ydtO 

ztO 

zdtO]; 

% 

%  initial  range  information 
% 

rxO=<xtO-xmO); 

rx(l)=rxO; 

ryO=(ytO-ymO); 

ry(l>=ryO; 

rzOKz^zmO); 

iz(l>=rzO; 

nn0=sqrt(xni0^2  +  yni0^2  +  2m0''2); 
nn(l>siTnO; 


itO=sqrt(xtO'^2  +  yt0^2  +  zt(y^2); 
rt(l)=rtO; 

lO  s=s<pt((xt0'Xin0y^2  +  (ytO'ymO)^2  +  (zt0-zm0y^2); 
i(l)=iO; 

% 

%  initial  time 
% 

time(l)=0.0; 

% 

%  SIMULATE  THE  SYSTEM 

% 

for  (i=l  ±max-l) 

% 

%  calculate  missile  and  target  speeds 

% 

vm(i)=sqrt(ms(24y'2  +  ms(44y^2  +  nis(6,i)'^2); 
vt(i>=sqrt(ts(24y'2  +  ts(44y'2  +  ts(64y^2); 

% 

%  calculate  line-of-sight  angles 

% 

sigma_pitch(i)=atan2((ts(54>ms(5  4)),sqrt((ts(l  4)-tns(  1 4))''2.. 
-Kts(34>-ms(34)y'2)); 


signut_yaw(i)=atan2((ts(34)-ins(34)),(abs(ts(  1 4>ms(  1 4)));; 


sigma(:4)=[signuLpitch(i) 
sigma_yaw(i)  ]; 

% 

%  calculate  missile  and  targe:  line-of-sight  angles 
% 

sigmaiii4)itch(i)=atan2(ms(54),sqrt(ms(l4)^2  +  ms(34)^2)); 
ngmam_yaw(i)  =aun2(ms(34),ms(14)); 
sigmatjpitch(i)=atan2(ts(54)4«irt(ts(14V'2  +  ts(34y'2)); 
ligmaLyaw(i)  =atan2(ts(34),:5(14)); 


% 

%  update  seeker  head  states 
% 

beta(:,i+l)=iAis*beia(:4)  +  dels*sigma(:4); 

% 

%  set  Lp  seeker  head  angle  rate  vector 
% 

betad(;4)=[beta(24) 

beta(44)]; 

% 

%  calculate  missile  and  target  flight  path  angles 

gammam_pitch(i)=atan2(ms(64).sqrt(ins(24y'2  + 
gammain_yaw(i)  =j»tan2(ins(44)4ns(24)); 
gammat^xh(i)=atan2(ts(64)4tin(ts(24>'2 1-  ts(4  i/^2)); 
gafninacjaw(i)  ssatar2(ts(44).ts(24)); 

% 

%  update  autopilot  states 

% 

ganiinadm(;4't-l)sphu.*gamfrKlr>r  •'  drla*NR*be»ad(:4); 

% 

%  compute  musiJe  vclof  .’v;  in  d***  RY  (pitcn,i  plane  (R’=Rxy) 

% 

if  (sigKUt4ntch(i>'=C.O; 

vm_pitch(i>ari»s(yin(i>^':cs(*igm  >  •‘W  '■.)-gidtunan\_yaw(i))); 

% 

%  limit  pitch  accelenaon  to  20  g's 

% 

if  (abs(vm_pitch(i)*gammadin(l4))<=644.0) 
an\4>itch(i)ssvn\_pitch(i)*ganunadm(  1 4); 
else 

am_pitch(i)=:644.0*sign(vtit_pitch0)*gammadm(  1  .i)); 
end 


else 


IP  IP 


vinjritch(i)=0.0; 

ain_pitch(i>=0.0; 

okI 

% 

%  calculate  missile  pitch  acceleration  vector  components 
% 

amLpilch_xy(i>=abs(am_pitch(i)*sm(sigma_pitch(i))); 
xddip_pitch(i)  =>am_4Jitch_xy(i)*cos(sigma_yaw(i)); 
yddin_piich(i)  =-am_pitch_xy(i)*sm(sigma_yaw(i)); 
rddm_pitch(i)  =am_pitch(i)*cos(sigma_pitch(i)); 

% 

%  compute  knissile  velocity  in  the  XY  (yaw)  plane 
% 

if  ((sigii.am_yaw(i)  -  sigmni_yaw(i))-^.0) 

vm_yaw(i)=abs(vm(i)*cos(gammam_pitch(i))); 

% 

%  limit  yaw  accelersdon  to  approximately  20  g's 
* 

if(abs(vni_yaw(i)*gammadni(24))<=644.0) 

aiH_yaw(i)=vm_yaw(i)*gammadm(24); 

else 

am_yaw(i)=644.0*sign(vin_yaw(i)*gammadm(24)); 

end 

else 

vin_ytw(i)=0,0; 

am_yaw(i)sO.O; 

end 

* 

%  calculate  missile  yaw  acceleration  vector  components 
% 

xddm_ya’v<i)=-am_yaw(i)*sin(sigma_jaw(i)); 
yddm_yaw(i)3s  am_yaw(i)*cos(sigmn,.yaw(i)); 

compute  overall  missile  acceleradon  components 


jR 


% 

xddm(i)=xddm_4)itch(i)  +  xddm_yaw(i); 
yddm(i)=yddni_pitch(i)  +  yddm_yaw(i); 
zddni(i)=zddni_pitch(i); 

% 

%  compute  total  missile  acceleration  magnitude 
% 

am(i)=sqrt(xddm(iy^2  +  yddm(iy^2  +  zddm(i)'^2); 

% 

%  generate  missile  input  vector 

% 

um=[xddm(i) 

yddm(i) 

zddm(i)]; 

% 

%  up<i4te  missile  states 

% 

ms(:4+l)»,3him*ms(:4)  +  delm^um; 

% 

%  set  target  acceleration  components 

% 

if(r(i)<=0.0) 

xddt(i)=-63*32J2*sin(gammat_yaw(i)); 
yddt(i)=  6J5*32J2*cos(gammaoaw(i)); 
2ddt(i)=0.0; 

else 

xddt<i)=0.0; 
yddt(iH).0; 
zddt(i)=0.0; 
end 

compute  target  acceleration  magnitude 
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at(i)=sqrt(xddt(iy^2  +  yddt(i)^2  +  2ddt(iy'2); 

% 

%  generate  target  input  vector 
% 

ut=[xddt(i) 

yddt(i) 

zddt(i)]; 

% 

%  update  target  states 

% 

ts(:4+l)=phit*ts(:4)  +  deU*ut; 

% 

%  calculate  tune  to  go 

% 

vt_pitch(i)=abs(vt(i)*cos(sigma_yaw(i)-gainmat_yaw(i))); 
idot(i>=  vt_pitch(i)*cos(gainmia_pitch(i)-signia_pitch(i))... 

-  vm_pitch(i)*cos(gamniain_pitch(i>sigmaj)itch(i)); 
ttg(i)=i<i)/(-rdot(i)); 

* 

%  calculate  updated  range  information 

% 

r(i+l)=sqrt((ts<U+l>ms(14+l)y'2  +  (ts(34+l>-ms(34+l))^2. 
+  (ts(54+l)-ins(54+l)m 

nn(i+l)=sqit(ins(14+iy'2  +  ins(34+l)^2  +  ins(54+l)^2); 

rt(i+l)=sqrt(ts(l  4+l)^2  +  ts(34+l)^2  + 15(54+1)^2); 

rx(i+lHts(14+l)-ms(14+l)); 

ry(i+l)s<tS(34+l)  -  ms(34+l)); 

rz(i+r)=(ts(54+l)  -  ms(54+l)); 

% 

%  set  up  missile  and  target  trajectory  data  for  plotting  in  3-D 

missile(i,:>=[ms(14)  nis(34)  ms(54)l; 
targeted:)  =[ts(14)  ts(34)  ts(54)l; 

% 
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%  update  tune  vector 

% 

time(i+l)=time(i)  +  dt; 

% 

%  check  to  see  it  engagement  is  at  closest  point  of  approach 
% 

if(r(i)<r(i+ 1  )),bicak,end 

end; 

% 

%  PLOT  RESULTS 

% 

%  time  to  go 

% 

plot(r(l:i-l),ttg(l:i-l)),gri(LxlabelCRange  -  FEET), 
ylabeI(Time-to-go  -  SECONDS’);titlcCRANGE  vs  TIME  TO  GO’); 
pause,clg 
% 

%  range  data 

% 

plot(time4'),grid,xlabelCTIME'),ylabeI('F£ET) 
titieCRANGE  VS  TIME*) 
pause,clg 
% 

%  missile  and  target  velocity  information 

% 

plot(time(l  :i),vm),grid,xlabel(TIME'),ylabelCFEET/SEC) 

titleCMISSILE  VELOCITY  VS  TIME’); 

pause,clg 

plot(tinM!(:ia).vt),grid,xlabel(TIME’),ylabcl(TEET/SE(^ 
tiUcCTARGET  VELOCITY  VS  TIME’); 
pause,clg 
% 

%  missile  and  target  acceleration  information 

% 
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plot(tiine(l:i-l),am(ld-l)),grid^bel(TIMEO,ylabel(T^T/SEC^^^ 

titleCMISSILE  ACCELERATION  VS  TIME’); 

pause.clg 

plot(time(ld).at),grid^abel(T[ME’),ylabel(TEET/SEC''2*) 
titlcCTARGET  ACCELERATION  VS  TIME’); 
pause,clg 
% 

%  Eight  path  angles 
% 

|dot(time(l:i),gaminain_pitch),gri(Lxlabel(TIME’),ylabelCRADIANS’) 
title(’MSL  PITCH  FUGHT  PATH  ANGLE  VS  TIME’); 
pause,clg 

plot(tiine(  1  :i),gainmain_y  aw),grid^abel(TIME’),ylabel(’RADIAN  S’) 
titlcCMSL  YAW  FUGHT  PATH  ANGLE  VS  TIME); 
pause,clg 

plot(tifne(l  d-l),ganunadm(l  ,(1  :i- 
l))),grid^abel(TIME’),yhbelOlADIANS/SECO 
titlcCMSL  PITCH  FUGHT  PATH  ANGLE  RATE  VS  TIME’); 
pause,clg 

plot(tiine(l  d-l),gaininadin(2,(l  d- 
l))).grid^bel(TIME'),ylabelCRADIANS/SEC) 
titlcCMSL  YAW  FUGHT  PATH  ANGLE  RATE  VS  TIME’); 
pause^lg 
% 

%  seeker  head  angle  information 

plpt(tme(ld-l).beta(l/Ld-l)))^d^belCTIME’),ylabelCRADIANS’) 

tifleCSEEKER  HEAD  PITCH  ANGLE  VS  TIME’); 

pause^lg 

plot(tiine(ld4),beta(2.(ld-l))).grid^abelCnME’),yIibelCRADIANS/SEC) 

titleCSEEKER  HEAD  PITCH  ANGLE  RATE  VS  TIME); 

pause/;lg 

plot(tiine(l:i-l).bcta(3Xld-l))),grid,xlabelCTIME’),ylabelCRADIANS’) 
titleCSEEKER  HEAD  YAW  ANGLE  VS  TIME’); 
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pause, c!g 

plot(time(l  :i-l),beta(4,(l  :i-l))),grid,xlabel(TIME’).yiabelCRADL\NS/SEC) 

titleCSEEKER  HEAD  YAW  ANGLE  RATE  VS  TIME’); 

pause,clg 

end; 
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APPENDIX  B-FORWARD  TIME  AND  ADJOINT  PROGRAMS 
% 


%  Lukenbill,  F.  C.  20  November  1990 

%  This  program  simulates  the  forward  time  and  adjoint  models  with  a 
%  step  acceleration  evasion  maneuver. 

% 

% 

%  Initialize  Variables 
% 

tm=1.0; 

tsh^.l; 

a=lAm; 

b=2*(l/tsh); 

c=(lAsh)*(l/tsh); 

N=4.0; 

Vc=2500; 

Tf=6.0; 

ki=l/(Vc*Tf); 

HE=0.0*pi/180; 

Vm=2000.0; 

VmHE=Vm*sin(HE); 

Atgts2S7.6;; 

% 


%  Input  State  Matricies 
% 


0 

0 

0 


-c*ki 

0 


1 

0 

0 

0 

0 

0 


0  0 

1  0 

-a  0 

0  0 

0  -c 

0  0 


0 

0 


0 

0 


a*N* 


Vc  0 


1  0 

-b  c*ki 

0  0 


0 

0 

0 

0 

0 

1 


0 

0 

0 

0 

0 

0 


-VmHE 

6 

0 

0 

0 

0 
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0  0  0  0  0  0  0  Atgt  0 

0  00  0  00000 
0  0  0  0  0  0  0  0  0]; 

% 

B=[0 

0 

0 

0 

0 

0 

0 

1 

11; 

% 

C=[-l  0  0  0  0  1  0  0  01; 

% 

%  Input  Timing  Information 
% 

dt=0.01; 

kmax:=Tf/dt+l; 

% 

%  Set  up  impulse  function 
% 

impulseszeros(  1  Jonax); 
impulse(l>=l/dt; 

% 

%  Initialize  vectors 
% 

x=zerosf  Janax); 
y=zen>s(l^ax); 
timeszeros(  1  ,kmax); 

% 

%  Simulate  the  System  (Forward  lime  Solution) 

% 
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I 


for  (i=l  Jonax-l) 

ki=l/(Vc*vTf-time(i))): 

A=[0  1  0  0  0  0  0 

0  0  1  0  0  0  0 

0  0  -a  0  a*N*Vc  0  0 

0  0  0  0  1  0  0 

^♦Id  0  0  -c  -b  c*ki  0 

0  0  0  0  0  0  1 

0  0  0  0  0  0  0 

0  0  0  0  0  0  0 

0  0  0  0  0  0  0 

[PhifwdJ)elfwd]=c2d(A3.dt); 
x(:4+l)  =  Phifwd*x(;4)  +  De’fv'd*impulse(i); 
y(:4-rl)  =  C*x(:4+1); 

tiine(i+l  >=tiine(i)+dt; 
end; 

y(i) 

% 

%  Plot  Results 


0 

0 

0 

0 

0 

0 

Atgt 

0 

0 


-VniHE 

0 

0 

0 

0 

0 

0 

0 

01; 


% 

plot(tiine(l:i),y(l,(l:i)));title<'Fcfwatd  Time  -  8  g  Step  Acceleration*); 
xlabelCnme');,ylabelCy  miss’);grid;pause; 

%  Reinitialize  variabler^  and  vectors 
% 

dt=0.01; 

kmaxssTf/dt+l; 

impulse=zeros(14cmax); 

impulse(l)=l/dt; 

x=2eros(94anax); 

y=2eros(14anax); 

tiineszeros(  1 4cniax); 

time(l)=dt; 

% 


%  Simulate  the  System  (Adjoint  Solution) 
% 


(i=l:kmax-l) 

ki= 

rl/(Vc*(time(i))); 

A^[0 

1 

0 

0 

0 

0 

0 

0 

-VmHE 

0 

0 

1 

A 

0 

0 

0 

0 

0 

0 

0 

0 

-a 

0 

a*N*Vc 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

-c*]d 

0 

0 

-c 

-b 

c*ki 

0 

0 

0 

0 

0 

0 

0 

0 

C 

1 

c 

0 

0 

0 

0 

0 

0 

0 

0 

Atgt 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0  0  0  0  0 
[Phiadj,Deladjl=c2d(A',C,dt); 

0 

0 

0 

0]; 

x(:,i+l)  *  Phiadj*x(:,i)  Deladj*impu!se(i); 
y(:4+l)  *  B’*x(:,i+1); 
time(i+l  )=time(i)+dt; 
end; 

Ai) 

%  Plot  Results 
% 

plot(time(l:i),y(l,(l:i)));title(' Adjoint  -  8  g  Srtp  Acceleration'); 

xl^!CTime');,ylabel('y  miss');gild;pause; 

end; 


99 


% 

%  Lukenbill,  F.  C.  20  November  1990 

%  This  program  simulates  the  forward  time  and  adjoint  models  with  a 
%  Barrel  Roll  maneuver. 

% 

% 

%  Initialize  Variables 


O', 


tm-1.0; 

tsh=0.1; 

a=:l/tm; 

b=2*(l/tsh); 

c=(l/tsh)*(l/tsh); 

N=4.0; 

Vc=2500; 

Tf=6; 

ki=V<Vc*Tf); 

HEs=u.O*pi/180; 

Vin=2000.0; 

VmHE=tVm*sin(HE); 

Atgt=:193.2; 

wt=1.0; 


^  =swt*wt; 


ATGT=wt*(Atgt*Atgt)/Tf; 


%  Input  State  Matricies 
% 


A={  0  10 

0  0  1 

0  0  -a 

0  0  0 

-c*ki  0  0 

0  0  0 

0  0  0 


0  0  0  0 

0  0  0  0 

0  a*N*Vc  0  0 

0  10  0 
-c  A,  c*ki  0 

0  0  0  1 

0  0  0  0 


0  -VmHEO  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  1  0 
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0  0  0  0  0  0  0  0  0  0  0 
0  ooooooooco 
0  0000000  00  1 
0  0  0  0  0  0  0  ATGT  0  -W  0] 

% 

B=[0 

0 

0 

0 

0 

0 

0 

1 

1 

0 

01; 

% 

C=[-l  dOOOlOOOOOl; 

% 

%  Input luning  Infoimadpn 
% 

dt=0.01; 

kmax=Tf/dHl; 

% 

%  Set  up  intpulse  function 
% 

impulse=zeros(  1  Jonax); 
iinpulse(l)=l/dt; 

%  , 

%  Initialize  vectors 
% 

xszeros(l  ljunax): 
yszeros(ljcinax); 
timeszeros(l  Janax); 
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% 

%  Stmulats  the  System  (Forward  Time  Solution) 

% 

for(i=ldanax-l) 

ki=l/(Vc*(Tf-time(i))); 

A=:(  0  1  0  0  0  0  0  0  -VmHE  0  0 

0  0  10  0  0  0  0  0  0  0 
0  0  -a  0  a*N*Vc  0  0  0  0  0  0 

0  0  0  0  1  0  00  0  0  0 

0  0  -c  'b  c*ki  0  0  0  0  0 

0  00  000  1  0000 
0  0  0  0  0  0  00010 
0  000  0  0  00000 
0  0  0  0  0  0  0  00  0  0 

0  co  o  0  0  0  0  0  0  1 

0  0  0  0  0  0  0  ATGT  0  -W  01; 

[PhifvrrtJDelfwd]=c2d(A,B40*. 
x(:d+l) »  Phifwd*x(:4)  +  Delfwd*impulsc(i); 
y(:,i+l)  =  C*x(:4+l); 
tiine(i+l  )=tiine(i>Hlt; 
end; 

y(i) 

* 

%  Plot  Results 
% 

plot(tiiTie(l:i).y(Uln)));tiUe(Ton»rard Time  -  6  g  Bane!  Roll'); 
xlabcI(Time*);,ylabelCy  miss’);grid;pausc; 

% 

%  Reinitialize  variables  and  vectors 
% 

dtsO.Ol; 
kmaxsTf/dt-fl; 
impulseszeros(  1  Janax); 
ijnpulse<l)sl/dt; 
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x=zcros(llJaTUx); 
y=zeros(lJanax); 
tiine=:zeros(  1  Jonax); 
tiine(l)=dt; 

% 


%  Simt'late  the  System  (Adjoint  Solution) 
% 

for  (i=ldunax-l) 


ki= 

:1/(VC< 

•(time(i))); 

> 

II 

O 

1 

0 

0 

0 

0 

0 

0  - 

VmHE 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-a 

0 

a*N*Vc 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

-c*ki 

0 

0 

-c' 

-b 

c*ki 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

ATGT 

0 

-W 

01 

[PhiadjJ>eladj]=<2d(A'.C.dt); 
x(;4+l)  =  Phiadj*x(:4)  +  Deladj*impulse(i); 
y<:4+l)  =  B**x(:4+l); 
tiiTie(i-fl  )=tunc(i)4dt; 
end; 

y<i) 

% 

%  Plot  Results 
% 

plot(time(l:i).y(].(ld)));titleC Adjoint  •  6  g  Barrel  Roll*): 

xlabclCTinfw’);,ylabelCy  miss*);grid;pause; 

end; 


103 


* 

%  Lukenbill,  F.  C.  21  November  1990 

%  This  program  simulates  the  forward  time  and  adjoint  models  with  a 
%  Spilt  'S’  maneuver 
% 

%  Initialize  Variables 

tm=1.0; 

tsh=0.1; 

a^tm; 

b=2*(l/tsh); 

c=(l/tsh)*(I/tsh); 

N=4.0; 

Vc=2500; 

Tf=<; 

ki*l/(Vc*T0; 

HE=0.0*pi/180; 

Vm=2000.0; 

VmHE=Vm*sin(KE); 

Atgt*193^; 

L=1.0; 

al=pi''2yL>'2; 

a2=pi*16*Atgt^2/(L*pi''2*T0; 

a3=9*al; 

i4sa2; 

% 

%  Initialize  State  Matricies 
% 

As[  6  l  0  0  0  0  0  0  -VmHE  0  000 

0  01  0  0  0  0  000  000 
0  0  -a  0  a*N*Vc  0  0  0  0  0  000 

0  0  0  0  10  0  00  0  000 
-c*ki  0  0  -c  -b  c*ki  0  0  0  0  000 

0  0  0  0  0  0  1  0  0  0  UOO 
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0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

010 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

100 

0 

0 

0 

0 

0 

0 

0 

^2 

0 

•al 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

001 

0 

0 

0 

0 

0 

0 

0 

a4 

0 

0 

0  -a3  01; 

% 

B=[0 

0 

0 

0 

0 

0 

0 

1 

1 

0 

0 

0 

01; 

% 

C=[-l  0  0  0  0  1  0  0  0  0  0  0  01; 
% 

%  Input  Uniing  Information 
% 

dt=0.01; 

kmaxsTf/dt-t-I; 

% 

%  Set  up  impulse  function 
% 

impulse=zeros(  1  Jcmax); 
impulsc(l>:l/dt; 

% 
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%  Initialize  vectors 
% 

x=zeros(  1  SJonax); 
y=2eros(14anax); 
ti2ne=:zeros(  1  Janax); 

% 

%  Simulate  the  System  (Forward  Time  Solution) 
% 

for  (i=l±max-l) 


ki=l/(Vc*(Tf-time(i))); 


> 

il 

o 

1 

0 

0 

0 

0 

0 

0 

-VmHE  0 

000 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

-a 

0 

a*N*Vc 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

000 

-c*ld 

0 

0 

-c 

-b 

c*ki 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

010 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

100 

0 

0 

0 

0 

0 

0 

0 

a2 

0 

-al 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

001 

0 

0 

0 

0 

0 

0 

0 

t4 

0 

0 

0  -a3  01; 

[Phifw(U>elfwd]ac2d(A3.dt); 
x(:4+l)  =  Phifwd*x(:4)  +  De!fwd*impalse<i); 
y(:,i+l)*C*x(:4+l); 
timefi+l  )=time(i>Kit; 
end; 
y(i) 

% 

%  Plot  Results 
% 

plot(time(l:i),y(l,(ln)));titleCForward  Time  ^  6  g  Split  S’); 
xlabel(Tlme*);,ylabel(’y  miss’);grid;pause; 
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% 

%  Reinitialize  variables  and  vectors 
% 

dfc=0.01; 
kinax=:Tf/dt-(-l; 
impulses:zeros(  1  dtmax); 
impulse(l)=l/dt; 
x=2eros(13Janax); 
y=2cros(l,kinax); 
time=zeros(  1  Jonax); 
time<l)=dt; 

% 

%  Simulate  the  System  (Adjoint  Solution) 
% 

for  (i=ldanax-l) 


ki=l/(Vc*(time(i))); 


> 

II 

o 

1 

0 

0 

0 

0 

0 

0 

-VmHE  0 

000 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

•a 

Q 

a*N*Vc 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

000 

-c*ki 

0 

0 

-c 

-b 

c*ki 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

.1 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

010 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

100 

0 

0 

0 

0 

0 

0 

0 

a2 

0 

•al 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

001 

0 

0 

0 

0 

0 

0 

0 

a4 

0 

0 

0-a3  0] 

(Phiadj,Deladj]=c2d(A*.C,dt); 
x(:4+l)  =  Phiklj*x(:4)  +  Deladj*impulse(i); 
y(:4+l)  =  B’*x(:d+l); 
tin^i+l  )=time(i)+dt; 
end; 


% 

%  Plot  Results 
% 

plot(time(ln),y(l,(l:i)));titleCAdjoint  -  6  g  SpUt  S’); 
xlabcl(Time');,ylabcl(’y  iniss');grid;pause; 
end; 


APPENDIX  C-DISSPLA  PLOTTING  PROGRAM 
CC  LukenbUl,  F.  C.  6/15/90 

CC  This  code  is  used  to  generate  the  three  dimensional  plots 
CC 

DIMENSION  MSLX(300).  MSLY(300).MSL4300), 
TG’DC(300),TGTY(300).TGTZ(300) 

NUM=100 
DO  10  I=1.NUM 

READ  (10,1000)  MSLXa).MSLYa).MSLZ(I) 

READ  (11,1000)  TGTXaXTGTYCD.TGTZd) 

10  CONTINUE 

1000  FORMAT  (3(1X,IE15.7)) 

CALLCOMPRS 
CALL  AREA2D(8.,8.5) 

CALL  X3NAME(’X  AXIS$M00) 

CALL  Y3NAMECYAXIS$M00) 

CALLZ3NAMECZ  AXISSMOO) 

CALL  VOLM3D(7.,7.,7.) 

CALL  VIEW(80000.,24000.,2000.) 

CALL  GRAF3D(0.,4000.,24000.,0.,4000,,24000.,0.,400.,2^ 
CALL  GRFrn(0.,0.,0.,0.,l  „0„0.,UL) 

CALL  AREA2D(7.,7.) 

CALL  GRi\F(0.,l.,6.,0.,l.,6.) 

CALL  GRID(2,2) 

CALL  END3GR(0) 

CALLGRFrn(0.,0.,0..0.,L,0.,L,L,0.) 

CALL  AREA2Cv;.,7.) 

CALL  GRAF(0.,4000.,24()00,0.,4000.,24000.) 

CALL  DASH 

CALL  CURVE  (MSLY>ISLXWM,0) 
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CALL  RESETCDASH) 

CALL  CURVE  (TGTY,TGTX.NUM.O) 

C  ALL  GRID(2,2) 

CALL  END3GR(0) 

CALL  DASH 

CALL  CURV3D(MSLXJ4SLY^SLX WM,0) 
CALL  RESETCDASH) 

CALL  CURV3D(TGTX,TGTY,TGTX,NUM,0) 
CALL  ENDPL(O) 

CALL  DONEPL 

STOP 

END 


no 


\ 
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